Epidemic thresholds and human mobility

A comprehensive view of disease epidemics demands a deep understanding of the complex interplay between human behaviour and infectious diseases. Here, we propose a flexible modelling framework that brings conclusions about the influence of human mobility and disease transmission on early epidemic growth, with applicability in outbreak preparedness. We use random matrix theory to compute an epidemic threshold, equivalent to the basic reproduction number \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_{0}$$\end{document}R0, for a SIR metapopulation model. The model includes both systematic and random features of human mobility. Variations in disease transmission rates, mobility modes (i.e. commuting and migration), and connectivity strengths determine the threshold value and whether or not a disease may potentially establish in the population, as well as the local incidence distribution.


A SIR metapopulation model
We consider a susceptible-infected-recovered (SIR, [15,3,8]) model describing the spread of an infectious disease over a discrete metapopulation, consisting of N patches connected by a weighted network. At each patch i, individuals remain at three possible states: susceptible, infected or recovered (denoted by S i , I i and R i respectively). These populations evolve according to the following system of differential equations where i ∈ {1, . . . , N } and N i = S i + I i + R i denotes the total population at patch i. Table 1 contains a description of the epidemiological parameters of the system. These parameters model the usual natural dynamics of human populations (births and deaths), as well as the life cycle of the disease (infection, recovery, loss of immunity). In addition to these, there are two qualitatively different terms in Eq.
(1) that model the mobility of individuals between patches.
Parameter Description Λ i Birth rate β i Transmission rate of the disease d i Natural death rate (excluding the modeled disease) δ i Disease-related death rate θ i Rate of loss of immunity to the disease α i Rate of disease overcome Table 1: Biological parameters of the model. Rates describing the natural dynamics of the population and the life cycle of the disease at patch i of the system (Eq. (1)).
The first of these terms corresponds to the short-term mobility between patches, and is labelled as commuting throughout the article. It models the infections generated at each patch i derived from interactions of susceptible individuals from patch i with infected individuals from other patches j. We assume that these new infections are generated at a rate β j c ji , that factors the transmission rate at patch j and a corrective factor c ji . Interactions may occur anywhere in the network and the corresponding individuals return to their patch of origin. We assume that the reference transmission rate is the one associated to the patch of origin j of the infected individual. The second term associated to human mobility concerns movement between patches with no return, and is labelled as migration. We assume that individuals go from patch j to patch i at a rate m ij and do not go back to their patch of origin. Depending on the chosen timescale for the system, this could represent daily or seasonal flows between locations in the network, for instance. Figure 1: The model. Three different processes may modify the number of infected individuals at each patch i. First, the natural dynamics of the disease, consisting of contagions of susceptible from infected individuals from the patch (at a rate β i of transmission of the disease) and depletion of infected individuals (at a rate γ i , that combines the rate of true recovery from the disease and the natural death rate of the population). Second, contagions of susceptible from infected individuals from other patches j (at a rate β j c ji , that factors the local disease transmission rate β j and a corrective term c ji quantifying the force of infection on patch i generated by interactions with individuals from patch j). Third, arrivals/departures of infected individuals from/to other patches j (at rates m ij and m ji respectively).
As mentioned in the text, both our commuting and migration terms are standard in metapopulation epidemic models [14,23], and are included to display the potential of RMT techniques in this context. Our choice of model is particularly suited for including spatial heterogeneities in transmission, such as disease strains, vaccination, or non-pharmaceutical interventions, but our methods could be developed analogously for other variations on the modelling of the transmission events (see for instance [18,22,17,4,30]).
Rather than choosing specific values for the coefficients c ij and m ij in Eq. (1), we assume that these are drawn from probability distributions. More precisely, we assume that every coefficient c ij is a realization of an iid random variable following a certain positive distribution c, with given mean and variance (µ c , σ 2 c ). Analogously, we assume that the coefficients m ij are drawn from iid random variables following another given positive distribution m, with mean and variance (µ m , σ 2 m ). As we show in Section 1.1 below, these two parameters (mean and variance of the distributions) are the only ones relevant for the possibility of spread of the disease over the system. No further constraints are assumed for the distributions, besides their positiveness. From a theoretical point of view, this flexibility should allow us to capture real-world mobility networks in our model, by simply computing the mean and variance of the empirically observed connectivity flows between the locations considered in the system. We incorporate nevertheless in Section 3 some structural variations in the network and show how these may affect the spread of the disease.
We will focus on the computation of epidemic thresholds, that determine the initial growth of the disease. These depend on the transmission and recovery rates of the disease and the mobility features of the network. Numerical simulations show that the remaining epidemiological and biological parameters in Eq. (1) influence instead the long-term behaviour of the system (Fig. 2). For instance, even an unstable scenario in which the infected population grows initially may result in a non-endemic disease if the disease-related death rate is high enough. In general, increasing the mortality rates produces an instant decay in the number of infected individuals, while if the natural death rate is lower than the birth rate and the system is unstable then the number of infected individuals will grow exponentially. An increase on θ, the rate of loss of immunity, results in a larger number of infected individuals at the endemic state of the disease. On the contrary, rising the rate of disease overcome α causes an endemic state with a smaller number of infected individuals. In both cases, the qualitative behaviour of the outbreak remains the same (Fig. 2).
In the following, we assume that the system in Eq. (1) is closed, that is, that the total population remains constant over time. This is a standard assumption in epidemic models for diseases with a relatively low death rate during their initial growth period [10]. This condition is achieved by assuming equal birth and natural death rates across patches (Λ i = d i for all i), and no diseaserelated mortality (δ i = 0 for all i). We denote the resulting recovery rate as γ i = α i + d i here and in the main article, combining the true recovery rate and the natural death rate. Under these hypothesis, the domain R 3N + is an invariant subspace under the action of the dynamical system (Eq. (1)), and thus the model is well-defined. This means in particular that if the initial conditions of the system are biologically feasible, then the system will remain to be such throughout time (see for instance Proposition 2 in [25] for a similar model).
Another consequence of the closedness of the system is that the disease-free state becomes in fact an equilibrium point of the dynamical system. The conditions I i = R i = 0 for all i characterize the disease-free equilibrium (DFE) of the system Eq. (1), that takes the form (S * 1 , 0, 0, S * 2 , 0, 0, . . . ). Its explicit expression can be found by solving the system of equations (see Theorem 3.3 That is, the distribution of the migration flows determines the distribution of the population over the patches of the system at the DFE, as the outgoing and incoming flows need to compensate in order for each local subpopulation to be constant. Note that the system (1) models the simplifying assumption that migration rates are independent of the state of the individuals.
The rightmost eigenvalue of the Jacobian matrix of the system at the DFE, s(J), governs the possibility of spread of the disease [9,10]. As follows from classical dynamical systems theory, if there is any eigenvalue with real part greater than zero the system is unstable at the DFE, i.e., the number of infected individuals tends to grow initially. If all the eigenvalues of the Jacobian at the DFE have negative real part, the system is stable and the disease dies out upon the introduction of a small number of infected individuals in the network. This behaviour provides the threshold condition analogous to that of R 0 exploited over the text. Moreover, the magnitude of the total 0.0e+00 5.0e+04  number of infected individuals over the system, the maximum of infected individuals at a single patch, and the time to arrive to the epidemic peak also depend non-linearly on the epidemic threshold s(J) (see Fig. 3 for some examples). Let us note that, unlike for systems modelling ecological communities [2,5], the differential equations for the susceptible, infected and recovered populations in (1) are coupled; it is only due to the known properties of SIR-like models [10] that the subsequent analysis by means of the Jacobian of the simpler infected subsystem is possible.

Random matrices
We assume for the moment that the rates of disease transmission and recovery of infected individuals are the same at all the nodes of the network, and denote them by β and γ respectively. Under these conditions, the Jacobian matrix of the infected subsystem of the model (1) at DFE reads as follows Recall that the commuting and migration coefficients are drawn from random variables with mean and standard deviation (µ c , σ c ) and (µ m , σ m ), respectively. As we argue below, the asymptotic eigenvalue distribution of J coincides with that of where 1 N is a N × N matrix of ones, the N × N matrix G N (0, σ) has iid entries following a certain distribution with mean zero and standard deviation σ = β 2 σ 2 c + 2βτ + σ 2 m (where τ = σ c σ m Cor(c ij , m ij )), and I N denotes the identity matrix of size N . We have chosen the monikers in Eq. (4) to help identify the role of each matrix in the eigenvalue distribution of J ′ : the structure matrix provides the mean of the off-diagonal entries of J ′ , the noise matrix captures their fluctuations around this mean, and the multiple of identity sets the mean of the diagonal.
We use the low rank perturbation theorem of Tao [26,21], that generalizes the well-known circular law from RMT [27], providing the asymptotic eigenvalue distribution of matrices such as J ′ whenever the structure matrix has low rank. Using this result (we actually use theorem 2.8 in [21], as we explain below), we find that the eigenvalue distribution of J ′ consists on two regions in the asymptotic N → ∞ regime. These are a circle on the complex plane on which the majority of the eigenvalues are uniformly distributed (the bulk), and a single outlier eigenvalue located on the real axis. After rescaling these distributions in N to match our model (see below for technical details), we find that the circle has radius σ √ N and is centered at β(1 − µ c ) − γ − N µ m , and that the outlier eigenvalue is located at βµ c (N − 1) + β − γ.
We now provide technical details that justify why the eigenvalue distribution of the original Jacobian matrix J given in Eq. (3) coincides in the asymptotic regime with the one of J ′ : • A central phenomenon in RMT is the universality of the asymptotic eigenvalue distributions of random matrices [27,12]: eigenvalues of matrices belonging to a large class of distributions share a common macroscopic behaviour in the large-N limit, that only depends on the mean and variance of the entries of the random matrix and not on the particular distribution followed by these. This is true for the circular law and also for the generalizations of it used throughout the article. As a consequence, the limiting distribution of the eigenvalues of J coincides with that of matrix J ′ for any choice of noise matrix, as long as its entries have the same variance σ 2 . Said differently, the eigenvalues of J follow the same asymptotic distribution for any choice of mobility coefficients having the same mean and variance and regardless of their particular values.
• Rigorous statements of the circular law [27], the low-rank perturbation theorem [26] and theorem 2.8 in [21] involve a different scaling of the random matrix J ′ , providing a fixed limiting distribution for the bulk of the eigenvalues. The matrix J ′ chosen here results in a rescaling of the distribution given by the low-rank perturbation theorem, after the substitutions µ → √ N µ, σ → √ N σ (for example in theorem 2.8 in [21]). The obtained distributions are also accurate in this case, as has been observed in previous works [2,1,13,24]. One of the advantages of RMT results is that often asymptotic convergence is apparent already at small values of N , and in particular we observe that is the case after this rescaling. See Fig. 4 for some examples of this fact, that holds for J as well.  The off-diagonal entries of the Jacobian have mean βµ c + µ m = 0.045 and standard deviation σ = 0.01 and the transmission and recovery rates verify β − γ = −0.4. Note that the predicted distribution is accurate for the outlier as well, even after the re-scaling in N explained in the text.
• Another difference between J and J ′ is that their diagonal elements are only equal in mean. O'Rourke and Renfrew prove [21] that the limiting eigenvalue distribution provided by the circular law is independent of the variance of the diagonal entries of the matrix, as long as this is finite. Note, however, that the results of RMT hold asymptotically for matrices with fixed parameters, and the variance of the diagonal entries of J grows with N . Numerical experiments show nevertheless a very close resemblance between the empirically computed eigenvalues of J and their theoretical expected distribution whenever σ m is relatively small. Even for larger values of this parameter, that controls the variance (N − 1)σ 2 m of the diagonal entries of J, the location of the outlier eigenvalue seems to be accurately predicted by the low-rank perturbation theorem, while the circular distribution of the bulk of the eigenvalues is deformed (see Fig. 5). For simplicity, we consider throughout the article parameter values for the random matrices such that the variance of their diagonal elements does not cause a significant perturbation of the circular law; these seem to match realistic assumptions on human mobility.
As explained in section 1, the system will be stable whenever all the eigenvalues of the Jacobian matrix J lie on the left half of the complex plane, and therefore one should in principle control for both the bulk region and the outlier eigenvalue to verify this condition. This results in the following two conditions for the stability of the system However, we see that if the first condition holds but the second is violated, then or, equivalently, that the coefficient of variation of the off-diagonal entries of J verifies Note, however, that the the off-diagonal entries of J follow a positive distribution, as explained in section 1. Without loss of generality, if we assume this distribution to be compactly supported on an interval [0, R] for some R > 0, the Bhatia-Davies inequality [7] introduces another constraint on the size of its variance where we have denoted by µ = βµ c + µ m the mean of the off-diagonal entries of J. Combining Eqs. (8) and (9), we see that in order for the first stability condition to hold and for the second to be violated, the distribution of the off-diagonal entries of J would need to verify In the asymptotic regime, this results in a contradiction, and even for finite but moderately large values of N this results in highly degenerated distributions that do not seem to model realistic systems. This means that for practical purposes, the outlier eigenvalue always lies on the right of the circle where the bulk of the eigenvalues is contained. Therefore, the stability of the system is determined by a single condition (the location of the outlier eigenvalue), which provides the desired epidemic threshold (Eq. (2) in main text)

Next-generation matrix
The basic reproduction number R 0 defines the number of contagions caused by a typical infected individual in a fully susceptible population. In compartmental models, it is given by the spectral radius of the next-generation matrix [9,28]. In this section, we show that RMT could be used to derive R 0 and s(J) for a simplified version of Eq. (1). Both numbers give equivalent epidemic thresholds; this is proved in [10] from the following equality The equivalence between the two epidemic thresholds follows from Eq. (12): s(J) < 0 ⇐⇒ R 0 < 1 : the disease will die out eventually (stability), s(J) > 0 ⇐⇒ R 0 > 1 : the disease may spread over the system (instability).
Comparing both conditions, we see that R 0 < 1 holds if all the eigenvalues of the next-generation matrix are contained in the unit circle ( Fig. 6), while s(J) < 0 if the maximum real part of the eigenvalues is smaller than zero.
We consider a simplified model with m ij = 0 in Eq. (1) and same transmission and recovery rates at the patches, given by and show that we can compute s(J) and R 0 with RMT, as well as its equivalence as presented in (13). In order to compute R 0 , we follow the algorithm presented in [10] for compartmental models, where R 0 is given by the spectral radius, ρ(K) of the next-generation matrix K. The next-generation matrix is computed as K = −T Σ −1 , where T and Σ can be extracted from the Jacobian matrix as J = T + Σ. Here, T is the transmission matrix, containing the terms related to production of new infected individuals, and Σ is the transition matrix, containing information relative to changes between infectious states. This decomposition must be chosen in such a way that T is a positive matrix and Σ is a positive off-diagonal matrix with negative dominant eigenvalue. For the infected subsystem (14) we can write and Σ = diag(−γ, ..., −γ). The next-generation matrix is then given by Following the same reasoning as in section 1.1, we find from RMT results that the eigenvalues of K are uniformly distributed in a circle with radius (β/γ)σ c √ N and center (β/γ)(1 − µ c ), with an outlier produced by the low-rank matrix µ c 1 N given by (β/γ)(µ c (N − 1) + 1). Therefore, constraining all these eigenvalues to lie inside the unit circle, we find that the basic reproduction number is given by If we compute the epidemic threshold s(J) from the Jacobian matrix, we obtain the following Equations (16) and (18) present R 0 and s(J), computed using the low-rank perturbation theorem [21] from RMT. Replacing these expressions in the inequalities (13) for the epidemic threshold and the basic reproduction number we see that An analogous reasoning holds for the case when R 0 > 0, showing the equivalence between both epidemic thresholds s(J) and R 0 .

Gershgorin and Perron-Frobenius theorems
We outline in this section how the classical Gershgorin [29] and Perron-Frobenius [19] theorems may be exploited in the analysis of our model. The less precise estimates provided by these results motivate further our use of RMT techniques. Let us first consider a model with no short-term mobility between patches. That is, we set all the commuting coefficients c ij to zero in (1), so that the infected subsystem simplifies to Around the DFE, the Jacobian of this system of equations reads Eigenvalue distribution for the Jacobian J and the next-generation matrix K of a stable system (left). The red dashed lines show the stability boundaries for the epidemic threshold s(J) (vertical line) and the reproduction number R 0 (unit circle). Eigenvalue distribution for an unstable system (right). The Jacobian has one eigenvalue with real part greater than zero. Accordingly, one eigenvalue of the next-generation matrix of the system is larger than 1 in modulus.

13
Using Gershgorin's circle theorem on the columns of this matrix, we see that its eigenvalues are contained in the union of the disks for i ∈ {1, . . . , N }. Therefore, a sufficient condition for the stability of the system is that the transmission rates β i are smaller than the exit rates γ i at all patches i, regardless of the values taken by the mobility coefficients m ji . Indeed, this is a consequence of the fact that infected individuals only move between patches in this case, without generating new infections in the process. Therefore, the only quantities that influence the growth of the disease are the rates β i and γ i at which new infections are being generated and depleted at each patch. An analogous phenomenon has been observed in [25] in the context of plant epidemics, see also [16] for a different proof of this fact. Thus, the stability of the system is determined by the local growth rates β i − γ i at the nodes of the network. We can also consider a system with nonzero commuting flows and bound the dominant eigenvalue of the Jacobian matrix using Perron-Frobenius theorem. For the simplest case of constant transmission and depletion rates across nodes, the theorem provides the bounds in terms of the entries J ij of the Jacobian matrix given in Eq. (3). These locate the epidemic threshold in an interval containing the point β − γ + βµ c (N − 1), with the width of the interval depending on the variance of the commuting rates σ c . Already in this case the RMT techniques provide a better estimate, providing an exact expression for s(J). A more significant improvement is found for the several generalizations of the model described in the following sections below: while RMT techniques provide in all cases an exact, analytic expression for s(J), the bounds obtained from Gershgorin and Perron-Frobenius theorems become much less precise.

Transmission and disease distribution
In Section 1.1 we obtained the epidemic threshold for a system with the same transmission and recovery rates at all patches. In this section we relax this assumption and study more general scenarios. Let us assume that at k nodes of the network, the transmission rate is β * = β + α. This could model the effect of successful vaccination, efficient public health measures (α < 0) or massive gatherings (α > 0), for instance. Labelling these nodes as the first k ones, without loss of generality, the Jacobian of the system around the DFE equals Using the same reasoning as in section 1.1, the eigenvalue distribution of this matrix can be seen to coincide with that of matrix (4) upon replacing the structure matrix by where a total of k columns have been modified in (23). We assume that the number of perturbed nodes k is small compared to N , to ensure that the low-rank perturbation theorem [21] provides a factual estimate for the corresponding eigenvalue distribution. Now k outliers arise in the eigenvalue distribution of the Jacobian, stemming from the k nontrivial eigenvalues of the matrix (23), with the largest one being equal to (Eq. (9) in main text) The resulting expression for the epidemic threshold is more complicated than that for the system with equal transmission rates. We also find a dependence on the average migratory flow µ m , in addition to the expected dependence on the perturbation in transmission α and the number of affected nodes k (see Fig. 7 for examples). Interestingly, an increase in transmission in one patch to β * = β + kα produces a stronger variation in the stability than an equivalent increase in k patches to β * = β + α (Fig.8b). We also find that Eq. (24) gives a better prediction for the epidemic threshold than the one obtained from the circular law Eq. (11) when considering the average transmission rate of the perturbed system β ′ = (kβ * + (N − k)β)/N (Fig.8a).
In addition to this, we may also consider random fluctuations in the transmission rates across patches, drawing β i from a probability distribution with mean µ β and variance σ β . Numerical experiments show that the eigenvalue distribution for the system with random β i is remarkably close to that predicted with a fixed diagonal term, corresponding to the mean of the random vector, µ β . The epidemic threshold would thus read as follows This prediction seems to be less accurate as either the coefficient of variation or the variance of the distribution followed by the β i increases (Fig. 9). Larger values of N provide a better agreement, and the larger the variance of the off-diagonal terms of J is with respect to the variance of its diagonal terms, the better the predicted distribution matches the empirically found one. This could model the case where there are some patches with zero or close to zero transmission rate. Even in this case, the average µ β determines the stability of the system and the number of infected  Figure 8: Variations in transmission at several patches of the system. Accuracy of the low-rank perturbation theorem (left). A better estimate is obtained for the rightmost eigenvalue in comparison to that obtained with Eq. (11) by simply using the average transmission rate of the perturbed system (size N = 100). Comparison between the epidemic threshold after increasing the transmission rate in a single patch to β * = β + kα versus increasing it at k patches to β * = β + α (right), for β = 0.05 and α = 0.5. There is a greater change in the rightmost eigenvalue when the increase is localized at one node compared to k nodes.
individuals may grow at all the patches. This is consistent with the results in [11,25], where it is proved that if the network is fully connected then almost surely there is no equilibrium point having patches with no infected individuals.
In section 1.1, the recovery rates were also assumed to be constant. We can also relax this assumption and assume that the recovery rate γ i in each patch is drawn from a probability distribution with mean µ γ and variance σ γ . This results in the diagonal entries of the Jacobian matrix having finite variance, which is still under the assumptions of the low-rank perturbation theorem [21]. We therefore obtain the same epidemic threshold in Eq. (11) in this case, upon replacing γ by µ γ .
Let us note that two seemingly contradicting phenomena have been described in this section. First, a different epidemic threshold has been obtained for targeted perturbations in the transmission rates of the system at a small number of patches by means of the low-rank perturbation theorem. In particular, this expression caused a larger increase in s(J) whenever an equivalent rise in transmission was concentrated at fewer locations. Secondly, we have observed a very close agreement between systems with random variations in their transmission rates and those with equal rates across patches. These findings seem to be incompatible, as the concentrated and disseminated increases in transmission considered in the first case were equal in mean. This fact may be understood in terms of the different nature of the two perturbations of the system. The first change favours a more severe rupture of the otherwise statistically balanced distribution of the underlying random matrix, while the second generates smoother statistical fluctuations that are compensated more easily. This is a subtle distinction that will arise also when considering variations on the mobility networks.

Random mobility networks
The epidemic threshold (Eq. (11)) shows that the average commuting flow is the only parameter from the connectivity networks that may influence the possibility of spread of the disease on unperturbed systems. Other features of the mobility flows may cause qualitative changes in the evolution of the pandemic, relevant for adequate design of public health measures. As explained in the main article, networks with a more heterogeneous distribution of commuting and migration flows result in more heterogeneous epidemic landscapes, in which patches with both higher and lower incidences can be found, see Fig. 10.
An instance of a more heterogeneous setting results from considering sparse connectivity networks. The examples of system (Eq. (1)) considered until now displayed fully connected networks. Here, we explore the stability of the system with a sparse connectivity network, i.e. a proportion of the flows is set to zero. In order to do this, each of the coefficients in the connectivity matrix is multiplied by a Bernoulli random variable with mean p; this models the effect of removing a proportion p of the connectivity flows of the system. This scenario has the same epidemic threshold as a system where the mean of the migration and the commuting coefficients is multiplied by the probability of success of the Bernoulli variable, µ = p(βµ c + µ m ). The latter represents a uniform : Accuracy on predictions for random transmission rates β i . Relative error for the outlier eigenvalue predicted by Eq. (11) for a system of N = 200 nodes with random transmission rates. Top row shows relative error in terms of the standard deviation σ β and a fixed average µ β = 1 (top left) and in terms of the average µ β with a fixed standard deviation σ β (top right). Bottom row shows relative error for the system in terms of the standard deviation σ β (bottom left) and the average µ β (bottom right) of the transmission rates for a fixed coefficient of variation CV = σ β /µ β = 0.5. If we compare the different scenarios, we observe a linear increase in loss of accuracy whenever CV is constant for either σ β and µ β . On the other hand, the relative accuracy of the prediction when CV is not fixed seems to depend non-linearly on σ β and µ β (and thus on CV as well). We see that as the CV increases the accuracy of the prediction decreases. proportional reduction in the strength of the connectivity on all of the flows of the system. Since the outlier eigenvalues coincide for both systems, we find the same qualitative behaviour for the spread of the disease (Fig. 11). However, the connectivity matrices are distinct and thus the bulk of the eigenvalues changes its shape: the radius of the circle differs due to the increase in variance in the sparse network (Fig. 11). In particular, the larger variance of the sparse system results in a more heterogeneous distribution of the incidence across patches.

Restrictions in commuting flows
We have observed the influence of the average commuting flow µ c in the epidemic threshold (Eq. (11)) when this remains the same for all the coefficients in the network. Let us explore the influence in the epidemic threshold of structured modifications in the commuting network. For this, we follow the same reasoning as in section 2 for the case of variable transmission rates. More precisely, we replace the structure matrix in (4) by one that adequately models a given modified mobility network, and compute the corresponding outlier by means of the low-rank perturbation theorem [21]. We consider a perturbation µ c → µ c + ν on the average commuting flow on some flows of the network, which increases or decreases the strength of the connectivity depending on the sign of ν, and assume that a number k of patches or a proportionate number of connectivity flows are perturbed. As above, we assume k to be small compared to the size of the network N in order for the low-rank perturbation theorem [21] to hold. We first consider targeted scenarios, where there is a set of k patches with modified average commuting flows: these comprise all the outgoing (A), incoming (B) or incoming and outgoing (C) flows. In this case, rows and/or columns in the structure matrix are modified. We also Time Figure 11: Sparse connectivity networks. Mobility network (top), eigenvalue distribution (center) and evolution of the number of infected individuals in time (bottom) for a sparse system (left, orange) and a uniformly restricted scenario (right, blue), with N = 50 and a proportion p = 0.6 of flows set to zero. Accordingly, the fully connected system shows a 60% reduction in the connectivity flows. The outliers are in the same position, in both scenarios the number of infected individuals increases initially. The bulk of eigenvalues is different, with higher variance for the sparse matrix implying a larger radius for the bulk of the eigenvalues and thus a higher variability between nodes. study non-structured scenarios: varying randomly chosen coefficients in the connectivity matrix (D), introducing symmetry in this random choice by modifying flows in both directions (E), and reducing uniformly the average commuting flow throughout the network (F). These scenarios are modelled as follows: • Strategy A: outgoing flows at selected locations. The structure matrix in (4) is now where the first k columns of the matrix have been modified, and the last (N − k) remain unchanged. For this choice of structure matrix, the outlier of the Jacobian (3) is given by • Strategy B: incoming flows at selected locations. The structure matrix is precisely the transpose of the structure matrix for strategy A, and thus the corresponding epidemic threshold coincides with Eq. (27).
• Strategy C: incoming and outgoing flows at selected locations. The structure matrix in (4) is now where the first k rows and columns of the matrix have been modified and bottom-right (N − k) × (N − k) block remains unchanged. For this choice of structure matrix, the outlier of the Jacobian (3) is given by • Strategy D, E and F: randomly chosen flows. The structure matrix is now replaced by (β(µ c + νk/N ) + µ m )1 N . The location of the outlier is now provided by the same expression as in the unperturbed case (3), with the corresponding modified average commuting flow As explained above when considering sparse matrices, the average commuting flow for strategy D and F coincide if the number of perturbed flows is the same, and thus the resulting epidemic thresholds are also equal. Strategy E consists of modifying the mobility flows in both directions at half as many randomly chosen edges of the network. While this perturbation results in a non-zero correlation between the diagonally opposite elements of the random matrix, this is known to only modify the distribution of the bulk of the eigenvalues [21], and thus results in the same epidemic threshold as strategy D as well. See section 3.3 below for more complicated correlations in the mobility matrices and their effect on the spread of the disease. Fig. 12 shows variations in the epidemic threshold s(J) corresponding to strategies A-F above upon variations in several relevant parameters of the system.

Perturbed migration scenarios
To investigate the response of the system to perturbations in the migration network, we compare in Fig. 13 six scenarios of increased migratory flows analogous to those considered in the previous section. We consider increased incoming flows to particular patches as well as outgoing flows, and non-structured modifications of the network. These could represent scenarios of humanitarian crisis or seasonal movements between particular locations, for instance.
As expected, the possibility of spread of the disease remains the same upon these perturbations in the network (see Results in main article and section 1 above). However, sensible differences arise in the local distribution of infected individuals over the system. For Scenario A (B), that models increased incoming (outgoing) migratory flows at certain nodes, we find that a larger (smaller) number of infected individuals concentrate at these locations, as a consequence of the overall larger number of individuals moving to (out of) the patches. If both incoming and outgoing flows are increased (Scenario C), these compensate each other and no significant differences arise in the affected patches. Among the randomly perturbed scenarios, those introducing more variability in the network (D and F) result in more heterogeneous epidemic landscape across patches.

Generally correlated networks
Baron et al [6] provide an asymptotic expression for the eigenvalue distribution of a random matrix with general pair-to-pair correlations between their elements. These may modify the location of both the bulk and the outlier eigenvalues, and thus change the stability of the corresponding system. Let us review their result and its implications for our model.
Consider a random matrix which elements z ij are drawn from a probability distribution with mean µ/N and variance σ 2 /N . Following [6], we assume that elements in the random matrix may  Figure 13: Variations in the migration network. We consider a base scenario (size N = 40) of sustained disease spread (positive, close to zero value for s(J), in grey) and test six perturbed migratory scenarios analogous to those described for the commuting network in the main article, here for an increase in migratory flows from µ m = 0.0005 to µ * m = 0.0025. Migratory rates concerning a total of 6 nodes or an equivalent proportion of flows are considered in the figure. In each scenario, the top graph displays the strength of the migration flow from patch i to patch j in its (i, j)-th cell. The brightness of the color represents the strength of the flow, with white representing absence of movement. Each line in the bottom graph shows the evolution of the infected population at the patches, colored according to its average incoming and outgoing migratory flows. The three targeted scenarios (in shades of orange) consist of A: increased incoming migration, B: increased outgoing migration, C: increased incoming and outgoing migration at half of the nodes. The three random scenarios (in shades of blue) consist of D: increasing randomly chosen unidirectional flows, E: increasing half as many randomly chosen flows in both directions, F: uniformly increasing all the flows in the network. also display non-trivial correlations with other matrix coefficients located at specific positions. More precisely, matrix elements that share a common index verify 1 ρ = Cor(z ij , z ji ), Γ N = Cor(z ij , z ki ), Then, as N → ∞, the bulk of the eigenvalues of the matrix distribute uniformly on the ellipse centered at the origin with semi-axes (1 + ρ)σ and (1 − ρ)σ, following the well-known elliptic law [20]. There is also a single outlier eigenvalue, located at This coincides with the location given by the low rank perturbation theorem (first term in the sum), with a translation that depends on the value of the correlations ρ and Γ, as well as the mean and variance of the matrix coefficients z ij (second addend). Additionally, Baron et al [6] show that any pair-to-pair correlations between elements of the matrix other than those in (31) may be absorbed in the mean of the coefficients z ij in the large-N limit. Therefore, eigenvalue predictions for networks that display more complex inter-dependencies can be approximated with accuracy by the expression (32) as long as N is sufficiently large. Interestingly, we see from this expression that the correlations r and c in (31) do not modify the eigenvalue distribution of the random matrix.
The parameters describing the statistical variability of the Jacobian matrix in our model can be obtained from those of the commuting and migration coefficients. Indeed, the off-diagonal terms of J matrix can be expressed as a linear combination of two differently distributed random matrices Let us denote as follows the correlations between matrix elements in both of these matrices (note the different scaling in N with respect to (31))  34)), which concern outgoing and incoming flows at the patches (center and right, respectively). These correlations result in more homogeneous distributions across the columns and rows of the mobility matrix respectively, favouring the existence of "source" or "sink" locations. See Fig. 5 in main article for examples of networks displaying nonzero Γ correlation.
See Figure 14 for some examples of correlated networks and the meaning of each parameter in (31) in the context of mobility networks. The relevant parameters encoding the statistical variability of the joint mobility network result from combining those of the two random matrices in (33) σ = β 2 σ 2 c + 2βτ + σ 2 m , ρ = (ρ c β 2 σ 2 c + ρ m σ 2 m )/σ 2 , Γ = (Γ c β 2 σ 2 c + Γ m σ 2 m )/σ 2 , r = (r c β 2 σ 2 c + r m σ 2 m )/σ 2 , c = (c c β 2 σ 2 c + c m σ 2 m )/σ 2 , where τ = σ c σ m Cor(c ij , m ij ). We have assumed in (35) that the correlations between other pairs of elements coming from each random matrix in (33) are null for simplicity (for instance, Cor(c ij , m ji ) = 0); the corresponding generalization would follow analogously. Note that we have chosen a different scaling in N from that in Baron et al [6]. This produces a deformation of the ellipse on which the bulk of the eigenvalues is uniformly distributed, as shown in Fig. 5 of the main text. The location of the outlier eigenvalue is nevertheless correctly estimated after performing the corresponding rescaling in equation (32), following an analogous reasoning as the one in section 1.1 for the circular law and the low-rank perturbation theorem. After correcting for the different diagonal terms in our model, this results in the following expression for the epidemic threshold for systems displaying the correlations in Eq. (34) See Fig. 5 in the main article for examples of this outlier eigenvalue and the effect of correlations in the stability of the system. Whenever both correlations Γ c and Γ m are zero, the last term in Eq. (36) vanishes and one recovers the original expression for the stability threshold Eq. (11). In the more general case, the sign of Γ determines whether the correlations facilitate or prevent the spread of the disease over the network (Fig. 15). We see from (36) that the epidemic threshold also depends on the correlation ρ between mobility flows going in opposite directions, as well the average migration rate and the joint variability of the connectivity network (Fig. 16).  36)) in terms of the joint correlation Γ for several choices of transmission rate β and joint mean and variance of the mobility network (µ = βµ c + µ m , σ 2 = β 2 σ 2 w + 2βτ + σ 2 c ). Negative (positive) values for this shift mean that the epidemic threshold for the correlated system is smaller (larger) than that expected from the unperturbed expression Eq. (11). The relatively small change in stability caused by ρ can be seen from the fact that a change from ρ = −0.3 (solid lines) to ρ = 0.6 (dashed lines) produces little variation in the outlier eigenvalue. In all cases, the sign of Γ determines whether the correlated network will be more stable (positive shift in s(J)) or less stable (negative shift in s(J)). Changes in Eq. (36) upon variations in several parameters of the system: the correlation between incoming and outgoing flows at the patches Γ, the correlation between flows going in opposite directions ρ, the average commuting flow µ c and the joint standard deviation of the connectivity networks σ. Note the different y scales in the three graphs; parameters for the base system are N = 200, Γ = 0.295, ρ = 0.44, µ c = 0.12, σ = 0.0012. As observed for systems with only nontrivial ρ correlation between their connectivity flows, the effect of this parameter on the epidemic threshold vanishes asymptotically (Eq. (36)).